/* _NSPROJ_C_
 * 
 * SOLVES THREE DIMENSIONAL, INCOMPRESSIBLE NAVIER-STOKES
 * EQUATION WITH RECTANGULAR DOMAIN USING FINITE DIFFERENCING 
 * ON A STAGGERED GRID
 * UTILIZES THE PROJECTION METHOD.  THE CONSERVATIVE FORM
 * CONVECTION IS DISCRETIZED WITHOUT THE USE OF UPWIND TREATMENT; 
 * BOTH THE ADVECTION AND DIFUSSION TERMS ARE TREATED EXPLICITLY. 
 * THE CURRENT SETUP IS FOR A LID DRIVEN CAVITY PROBLEM WITH A WELL
 * KNOWN VELOCITY AT THE TOP AND NO-SLIP BOUNDARY CONDIITIONS AT
 * ALL OTHER BOUNDARIES
 * 
 * ORIGINALLY WRITTEN BY DR. ZHI-GANG FENG: UNIVERSITY OF NORTH TEXAS
 * 
 * PORTED TO C AND 3D BY JOHN LEIDEL: GRADUATE STUDENT
 * 
 */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <sys/types.h>
#include "nsproj.h"

//-----------------------------------------START_SIM
// START_SIM
//--------------------------------------------------
int start_sim( char *datafile )
{
    /* VARS */
    double Re          = 0;     //-- Reynolds Number
    double dt          = 0;     //-- Time Step
    double ct          = 0;     //-- Current Time
    double tf          = 0;     //-- Final Time
    double lx          = 0;     //-- Width of Box
    double ly          = 0;     //-- Height of Box
    double lz          = 0;     //-- Depth of Box
    double dx          = 0;     //-- delta_x
    double dy          = 0;     //-- delta_y
    double dz          = 0;     //-- delta_z
    double simtime     = 0;     //-- TOTAL SIMULATION WALLCLOCK TIME
    double lsimtime    = 0;     //-- TOTAL ITERATION WALLCLOCK TIME
    double step1       = 0;     //-- TOTAL TIME FOR STEP_1
    double step2       = 0;     //-- TOTAL TIME FOR STEP_2
    double step3       = 0;     //-- TOTAL TIME FOR STEP_3
    double iot         = 0;     //-- TOTAL TIME SPENT IN I/O
    double bt          = 0;     //-- TOTAL TIME SPENT IN BOUNDARY CONDITIONS
    int nsp            = 0;     //-- DETERMINES NSP FILE OUTPUT
    int tec            = 0;     //-- DETERMINES TEC FILE OUTPUT
    int nt             = 0;     //-- NUMBER OF ITERATIONS [TOTAL_TIME/dT]   
    int c              = 0;     //-- Counter
    

    struct timeval starttime;
    struct timeval endtime;
    NSMESH mesh;
    char project[1024];
    /* ---- */

#ifndef _OPENMP_
    printf( "%s\n", "NSPROJ [single threaded] v. 0.2" );
#else
    printf( "%s\n", "NSPROJ [multi-threaded: OpenMP] v. 0.2" );
#endif

#ifdef _GCC_
    set_fpu (0x27F);
#endif

    /* PARSE CONFIG FILE */
    if( nsproj_readdatafile( datafile, 
			     &Re, 
			     &dt, 
			     &tf, 
			     &lx, 
			     &ly, 
			     &lz, 
			     &mesh.uN,
			     &mesh.uS,
			     &mesh.vE,
			     &mesh.vW,
			     &mesh.wF,
			     &mesh.wB,
			     &mesh.mx, 
			     &mesh.my, 
			     &mesh.mz,
			     &nsp,
			     &tec,
			     project ) != 0 )
    {
	printf( "%s%s\n", "ERROR : COULD NOT READ DATA FILE : ", datafile );
	
	return -1;
    }

    /* SETUP GRID POINTS && TIMING */
    nt = (int)ceil( tf/dt);
    dt = tf/nt;
    dx = lx/mesh.mx;
    dy = ly/mesh.my;
    dz = lz/mesh.mz;

    printf("tf=%f dt=%f nt=%d\n",tf,dt,nt);

    /* SANITY CHECK VARIABLES */
    if( (dt<=0) || (dx<=0) || (dy<=0) || (dz<=0) || (lz<=0) || (ly<=0) || (lz<=0) || (nt<1) )
    {
	printf( "%s\n", "ERROR : DATA OUT OF BOUNDS" );
	
	return -1;
    }

    /* ALLOCATE MEMORY FOR MESH */
    if( nsproj_meshalloc( &mesh ) != 0 )
    {
	printf( "%s\n", "ERROR : COULD NOT ALLOCATE MEMORY FOR MESH" );
	
	return -1;
    }
    
    /* INIT MESH */
    if( nsproj_meshinit( &mesh ) != 0 )
    {
	printf( "%s\n", "ERROR : COULD NOT INITIALIZE MESH ELEMENTS" );
	nsproj_meshfree( &mesh );
	
	return -1;
    }

    printf( "%s\n", "....BEGINNING SIMULATION" );
#ifdef _OPENMP_
    /* TODO :: PRINT NUMBER OF OPENMP THREADS */
#endif

    printf( "%s\n", "....TIMING IN SECONDS" );
    printf( "\n%s\n", "ITER\t CUR_T\t\t STEP_1_T\t STEP_2_T\t STEP_3_T\t BC_T\t\t IO_T\t\t ITER_T\t\t TOTAL_T" );


    /* BEGIN SIMULATION */
    for( c=0; c<nt; c++ )
    {
	/* SOLVING TENTATIVE USTAR/VSTAR/WSTAR  WITH IMPLICIT METHOD */
	gettimeofday( &starttime, NULL );
	
	velstar( Re, dt, dx, dy, dz, &mesh );
		
	gettimeofday( &endtime, NULL );
#ifdef _DEBUG_
	printf("u=%f v=%f w=%f p=%f\n",mesh.U[3][3][3], mesh.V[3][3][3], mesh.W[3][3][3],mesh.P[3][3][3]);
#endif
	step1 += (double)(endtime.tv_sec-starttime.tv_sec);
	step1 += (double)(endtime.tv_usec-starttime.tv_usec)*0.000001;

	/* SOLVING PRESSURE EQUATION */
	gettimeofday( &starttime, NULL );
	
	pressure( dt, dx, dy, dz, &mesh );
		
	gettimeofday( &endtime, NULL );
#ifdef _DEBUG_
	printf("u=%f v=%f w=%f p=%f\n",mesh.U[3][3][3], mesh.V[3][3][3], mesh.W[3][3][3],mesh.P[3][3][3]);
#endif
	step2 += (double)(endtime.tv_sec-starttime.tv_sec);
	step2 += (double)(endtime.tv_usec-starttime.tv_usec)*0.000001;

	/* PRESSURE CORRECTION */
	gettimeofday( &starttime, NULL );
	
	pressure_correction( dt, dx, dy, dz, &mesh );
		
	gettimeofday( &endtime, NULL );

	step3 += (double)(endtime.tv_sec-starttime.tv_sec);
	step3 += (double)(endtime.tv_usec-starttime.tv_usec)*0.000001;


	/* BOUNDARY CONDITIONS */
	gettimeofday( &starttime, NULL );

	boundary_cond( dt, dx, dy, dz, &mesh );

	gettimeofday( &endtime, NULL );

	bt += (double)(endtime.tv_sec-starttime.tv_sec);
	bt += (double)(endtime.tv_usec-starttime.tv_usec)*0.000001;

	/* OUTPUT ALL CURRENT DATA */
	gettimeofday( &starttime, NULL );

	if( nsp )
	{
	    nsproj_writedata( project, mesh, ct );
	}
	if( tec )
	{
	    nsproj_writetecdata( project, mesh, c );
	}
	gettimeofday( &endtime, NULL );

	iot += (double)(endtime.tv_sec-starttime.tv_sec);
	iot += (double)(endtime.tv_usec-starttime.tv_usec)*0.000001;

	lsimtime = step1 + step2 + step3 + iot + bt;
	simtime += lsimtime;

	/* PRINT STATUS TO SCREEN */
	printf( "%d%s%f%s%f%s%f%s%f%s%f%s%f%s%f%s%f\n", c, "\t",ct, "\t",step1, "\t",step2, "\t",step3, "\t",bt, "\t",iot, "\t", lsimtime, "\t",simtime );

	/* RESET LOCAL TIMING */
	step1    = 0;
	step2    = 0;
	step3    = 0;
	lsimtime = 0;
	iot      = 0;
	bt       = 0;
	ct       = ct+dt; /* UPDATE CURRENT TIME [t] */
    }

    printf( "%s%s%s\n", "NSPROJ :: ", project, " COMPLETE!" );
    printf( "%s%f\n", "TOTAL SIMULTION WALLCLOCK TIME = ", simtime );

    
    nsproj_meshfree( &mesh );
    /* DONE */

    return 0;
}

//-----------------------------------------MAIN
// MAIN
//---------------------------------------------
int main( int argc, char **argv )
{
    /* VARS */
    char *datafile = NULL;
    int rtn        = 0;
    /* ---- */

    if( argc < 2 )
    {
	fprintf( stderr, "%s\n", "USAGE: nsproj <datafile>" );
	return -1;
    }

    datafile = (char *)malloc( 1024 );
    if( datafile == NULL )
    {
	fprintf( stderr, "%s\n", "ERROR :  COULD NOT ALLOCATE MEMORY FOR DATA FILE" );
	return -1;
    }

    sprintf( datafile, "%s", argv[1] );
    
    rtn = start_sim( datafile );
    
    free( datafile );
    datafile = NULL;

    return rtn;
}
